Proceedings Title IAU Symposium 

Proceedings IAU Symposium No. IAU Symposium No. 252, £@9<§008 International Astronomical Union 
A. C. Editor, B.D. Editor & C.E. Editor, eds. DOI: 00.0000/X000000000000000X 

Hydrodynamic simulations of the core 

helium flash 

Miroslav Mocak, Ewald Miiller 
Achim Weiss and Konstantinos Kifonidis 

Max-Planck-Institut fur Astrophysik, Postfach 1312, 85741 Garching, Germany 
email: mmocak@mpa-garching.mpg.de 

Abstract. We desribe and discuss hydrodynamic simulations of the core helium flash using an 
initial model of a 1.25 Mq star with a metallicity of 0.02 near at its peak. Past research concerned 
with the dynamics of the core helium flash is inconclusive. Its results range from a confirmation of 
the standard picture, where the star remains in hydrostatic equilibrium during the flash (Deupree 
1996), to a disruption or a significant mass loss of the star (Edwards 1969; Cole & Deupree 1980). 
However, the most recent multidimensional hydrodynamic study (Dearborn et al. 2006) suggests 
a quiescent behavior of the core helium flash and seems to rule out an explosive scenario. Here 
we present partial results of a new comprehensive study of the core helium flash, which seem 
to confirm this qualitative behavior and give a better insight into operation of the convection 
zone powered by helium burning during the flash. The hydrodynamic evolution is followed on 
a computational grid in spherical coordinates using our new version of the multi-dimensional 
hydrodynamic code HERAKLES, which is based on a direct Eulerian implementation of the 
piecewise parabolic method. 
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1. Introduction 

First results on the core helium flash were gained from one-dimensional hydrostatic 
numerical simulations of a 1.3 M star (Z=0.001) (e.g., Schwarzschild & Harm 1961). 
During the flash, the star underwent a thermal runaway due to the ignition of helium 
under degenerate conditions in its center. It reached a peak at maximum core temperature 
of ~ 3.510 8 if and total energy generation rate of ~ 10 12 L Q . The calculations were 
redone later with better numerical techniques and improved treatment of major physical 
processes (Sweigert & Gross 1978) and although the ignition of helium occured off- 
center due to neutrino processes, they did not change the general picture mentioned 
earlier. It turns out, that the typical e-folding times for the energy release from helium 
burning become as low as hours at the peak of the flash, and therefore arc comparable 
to convective turnover times. Thus, the usual assumptions used in simple descriptions 
of convection in one-dimensional hydrostatic calculations (e.g. instantaneous mixing) do 
not have to be valid any longer. Previous attempts to relax these assumptions by allowing 
for hydrodynamic flow remained inconclusive (Edwards 1969; Deupree 1996; Dearborn ct 
al. 2006). Using a modified version of the HERAKLES code (Kifonidis ct al. 2003) which 
is capable of solving the hydrodynamic equations coupled to nuclear burning and thermal 
transport in up to three spatial dimensions, we want to deepen our understanding of the 
convection during the core helium flash at its peak investigating it by means of two- 
dimensional and three-dimensional hydrodynamic simulations. 
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Figure 1. Left panel: Temperature (in 10 7 K, solid) and density (in 10 s g cm -3 , dotted) dis- 
tribution of the initial model M. Right panel: Chemical composition of the initial model M, 
showing dominant fraction of helium. 



Table 1. Some properties of the initial model: total mass M, stellar population, metal content 
Z, mass Mhe and radius Rue of the helium core (X( 4 He) > 0.98), nuclear energy production 
in the helium core Lhb, maximum temperature of the star T max , and radius r max and density 
Pmax at the temperature maximum. 
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2. Initial setup 

The initial model was obtained with the stellar evolution code GARSTEC (Weiss & 
Schlattl 2007) . Some of its properties are listed in Table[l] The temperature, density and 
composition distribution of the model is depicted in Figure[l] The model encompasses a 
white dwarf-like degenerate structure with an off-center temperature maximum resulting 
from plasma- and photo-neutrino cooling and a central density of about 710 5 g cm~ 3 . 
The isothermal region in the center of the helium core is followed by almost discontin- 
uous jump in temperature up to T max ~ 1.7 10 s K and convection zone driven by the 
superadiabatic temperature gradient. The model is composed mostly of helium 4 He with 
an abundance X( 4 He)> 0.98. The remaining composition of the stellar model is ^^H, 3 He, 
12 C, 13 C, 14 N, 15 N and 16 0. For our hydrodynamic simulations we adopt the abundances 
of 4 He, 12 C and 16 from the initial model, since the triple-a reaction dominates the 
nuclear energy production rate during the flash. The remaining composition is assumed 
to be adequately represented by a gas with a mean molecular weight equal to that of 
20 Ne. 



3. Hydrodynamic simulations 

Table[2] summarizes some characteristic parameters of our two-dimensional (2D) and 
three-dimensional (3D) simulations that are based on model M. They were performed 
on an equidistant spherical grid encompassing 95% of the helium core's mass except for 
a central region with a radius of r = 2 10 8 cm, which was excised in order to allow for 
larger timesteps. 
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Table 2. Some properties of the two and three-dimensional simulations: number of grid points 
in radial (N r ) and angular (Ng,N^) dimension, radial (Ar in 10 8 cm) and angular (AO, A0) 
resolution, characteristic length scale l c (in 10 8 cm) and velocity v c (in 10 6 cm s _1 ) of the flow, 
respectively, expansion velocity at the position of temperature maximum v exp (in cm s _1 ), 
entrainment rate v en t of the outer convective boundary (in m s _1 ), typical convective turnover 
time to and maximum evolution time t m ax (in s), respectively. 



run | N r x Ne x N^> Ar AO A<j> l c v c v exp v ent t a t max 

DV2 180 x 90 5.55 2.° - 4.7 1.03 - 6. 7. 910 30000 

DV4 360 x 240 2.77 0.75° - 4.7 1.52 +92. 14. 620 60000 

TR 180 x 60 x 60 5.55 1.5° 1.5° 4.7 0.7 +6. 7. 1340 5300 
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Figure 2. Snapshots of the evolved convection at 60000 s in model DV4, showing the velocity 
field (upper panel), and the velocity amplitude |v| in 10 6 cm s _1 (bottom panel), respectively. 

All our 2D and 3D models undergo initially (t < 1200 s) a common evolution where 
convection sets in after roughly 1000 s. During this phase, hot bubbles appear in the 
region where helium burns in a thin shell (r~ 510 8 cm). After ~ 200 s, they cover 
complete height of the convective region and reach a steady state with several upstreams 
(or plumes) of hot gas carrying the released nuclear energy away from the burning region, 
thereby inhibiting a thermonuclear runaway. 

Fully evolved convection (t > 1500 s) in 3D is significantly different than in 2D, since 
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Figure 3. Left upper panel: Convective energy flux Fc in model DV2 (dotted) and model TR 
(dash-dotted) and kinetic energy flux Fk in model DV2 (solid) and model TR (dashed) (Hurlburt 
et al. 1986). Right upper panel: The r.m.s convection velocity in model DV2 (dashed), DV4 
(dotted) and TR (solid) overplotted with the convective velocities predicted by mixing-length 
theory (long-dashed). Bottom left panel: Temporal evolution of the outer convective boundary 
in model DV2 (dash-dotted), DV4(dashed) and TR (solid), respectively. Bottom right panel: 
Expansion velocity v exp in model DV2 (dash-dotted), DV4 (dotted) and TR (solid) together 
with the expansion velocities of the initial stellar model (long-dashed). 



the shape of turbulent streams which transport energy is totally distinct. However, the 
amount of energy which needs to be transported by the convection in order to prevent 
a thermonuclear runaway during the flash is in both cases similar. The resulting typical 
convective velocities are therefore much higher in 2D than in 3D (Fig. [3]). 

The structural differences between 2D and 3D flows are clearly visible in the distribu- 
tion of the kinetic flux across the convection zone (Fig. [3]). The typical evolved 2D flows 
contain well defined vortices (Fig.[2]) with their central regions never interacting with 
the region of the dominant nuclear burning above the T max . This results in a reduced 
kinetic flux between 5 10 8 cm < r < 6 10 8 cm, since the gas in that region, on average, is 
located at bottom of convective vortices, and thus does not experience any strong radial 
flow. On the other hand, the distribution of the kinetic flux in 3D is rather smooth, and 
the flow structures tend to be also smaller than in 2D. This is apparent when comparing 
Figure[4] with Figure[2j The 2D structures (vortices) have an angular size of around 40°. 
The structures in the 3D are column shaped, with a smaller angular size. The convective 
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Figure 4. Left panel: Isosurface of a radial velocity field of model TR at t = 4150 s. The 
bright color marks a positive velocity of +6 10 5 cm s~ 1 (the upflow streams) and the dark 
color a negative velocity of — 6 10 5 cm s _1 (downflow streams). Axis tickmarks are in units 
10 8 cm. Right panel: Fractional volume occupied by upflow and downflow streams in model 
DV2 (dashed-dotted), DV4 (dashed) and TR (solid), respectively. 



and kinetic flux is lower in 3D than in the 2D, but the total energy production is about 
20 % higher in 3D (because no symmetry restrictions are imposed, and due to the strong 
dependence of the triple-a reaction rate on the temperature). The convective and kinetic 
flux carry together more than 90 % of the energy produced by the burning. The 3D 
velocities qualitatively match the velocities predicted by the mixing-length theory better 
than in 2D, where the velocities are clearly overestimated. Figure[3] shows that they also 
depend on resolution, being higher in the simulation with the highest resolution. 

The extent of the convection zone increases with time. Due to turbulent entrainmcnt 
(Meakin & Arnett 2007), convective boundaries defined by the Schwarzschild criterium 
are pushed towards the center of the star, and towards the stellar surface, respectively 
(Fig. [3]). This is in contradiction with the predictions made by (ID) hydrostatic stellar 
modeling. For ilustration, the temporal evolution of the location of the outer convective 
boundary is depicted in Figure[3] It is defined as the radius, where the mean carbon abun- 
dance X( 12 C ) ~ 0.002. The rapid initial jump of the boundary position to r ~ 9.4 10 8 cm 
at about ~ 1200 s is due to the first touch of the convective flow on the boundary. Later 
entrainment is rather steady. The velocity of the outer boundary (entrainment rate) in 
our models are listed in Table|2] The entrainment involves a few radial zones only over the 
longest simulation we have performed. Although the 12 C abundance distribution stayed 
discontinuous at boundaries (no evident effect of numerical diffusion is detected) , the en- 
trainment rates presented here have to be considered as an order of magnitude estimate 
only. The entrainment at the inner convective boundary occurs with a rate much smaller 
than at the outer convective boundary. Therefore it is not discussed further here, since 
longer simulations are needed for definite statements about its evolution. 

A similar feature which 2D and 3D seem to share is the upflow-downflow asymmetry. 
The downflows cover a much bigger volume in the convection zone than the upflows. The 
downflows dominate more in 3D. Interestingly, the kinetic flux is always positive in both 
cases, although the downflows fill almost the whole convection zone. This implies that 
the downflows are much slower then the upflows. 

The expansion velocities v exp = M r / '4irr 2 p are in good agreement with those of initial 
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Figure 5. Maps of 12 C abundance (in units of 10~ 3 ) in a horizontal plane in model TR at t — 
5231 s, at different radii: n = 4.8 10* cm (left), r 2 = 6.5 10 8 cm (middle), r 3 = 9.3 10 8 cm (right). 

stellar model only in the 2D model with highest resolution DV4 (Fig. [3]). The expansion in 
the low resolution models does not match at all. Due to the different dynamic properties 
of the flow in the less resolved models, the spherical mass flow is weaker. 

Mixing at the convective boundaries and across the convection zone in 2D and 3D is 
quite different as well. In 2D, due to the symmetry restriction, every turbulent feature 
is in fact an annulus. Hence, turbulence and mixing can be properly modelled only by 
means of 3D simulations. The most apparent turbulent structures during the flash, in 
3D, look at the bottom of the convection zone like thin hot fibers enriched by carbon 
and oxygen (ashes from the helium burning). The flow then gets more uniform inside 
the convective region, but looks more turbulent again at the outer convective boundary 
(Fig.§. 



4. Conclusions 

We find that the core helium flash neither rips the star apart, nor significantly alters 
its structure. The evolved convection in 3D looks different from that in 2D. Typical 
convective velocities arc higher in 2D than in 3D where they also tend to fit the predictions 
made by mixing length theory better. Hydrodynamic simulations show the presence of 
turbulent entrainment, which results in a growth of the convection zone on dynamic time 
scales. 
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